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procedure wofz(x,y,re,im) ; value x,y; real x,y,re,im; 

comment This procedure evaluates the real and imaginary part of the function 
2 

w(z) = exp(-z )erfc(-iz) for arguments z = x+iy in the first quadrant of 
the complex plane. The accuracy is 10 decimal places after the decimal 
point, or better. For the underlying analysis, see W. Gautschi, "Efficient 

fb ir cn 

computation of the complex error function," submitted fee SIAM J. Math. 

Anal. ; 

begin integer capn, nu, n, npl; real h, h2, lambda, rl, r2, s, si, s2, tl, 
t2, c; Boolean b; 
if y < 4.29 /\ x < 5.33 then 
begin 

s := (l-y/4.29) x sqrt(l~x x x/28.41); 
h := 1.6 x s; h2 := 2 x h; 


capn := 6+23 x s; nu := 9 + 21 x s 


end else begin h := 0; capn := 0; nu := 8 end ; 
if h > 0 then lambda := h2 + capn; 
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b := h = 0 V lambda = 0; 
rl := r2 ;= si := s2 ;= 0; 
for n ;= nu step -1 until 0 do 
begin 

npl :« n+1; 

tl °.= y+h+npl x rl; t2 :*» x-npl x r2; c ; = .5/(tl x tl + t2 x t2) ; 
rl := c x tl; r2 := c x t2; 
if h > 0 /\ n < capn then 
begin 

tl ; = lambda + si; si s= rl x tl - r2 x S 2; s2 := r2 x tl + rl x s 2 
lambda lambda/h2 
end 
end; 

re s= if_ y = 0 then exp(-x x x) else 

1. 12837916709551 x ( if b then rl else si); 
im 1.12837916709551 x (if b then r2 else s2) 


end wofz 
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Author : Walter Gautschi, Department of Computer Sciences, 

‘"Purdue University, Lafayette, Indiana 47907 

Abstract : The paper is concerned with the computation of w(z) = 

2 

exp(-z )erfc(-iz) for complex z «* x+iy in the first quadrant Q^: 
x >_ 0, y >_ 0o Using Stieltjes* theory of continued fractions it 
is first observed that the Laplace continued fraction for w(z) , 
although divergent on the real line, represents w(z) asymptotically 
for large z in the sector S: -ir/4 < arg z < 5ir/4„ Specifically, the 
n-th convergent approximates w(z) to within an error of 0(z as 

z <» in S» A recursive procedure is then developed which permits 
evaluating w(z) to a prescribed accuracy for any z e Q^„ The pro- 
cedure has the property that as |z| becomes sufficiently large, it 
automatically reduces to the evaluation of the Laplace continued 
fraction,, or, equivalently, to Gauss-Hermite quadrature of 


(i/ir) 


exp(-t )dt / (z-t) 


Key words and key phrases : error function for complex argument, Voigt 

function, Laplace continued fraction, Gauss-Hermite quadrature, re- 


cursive computation 



EFFICIENT COMPUTATION OF THE COMPLEX ERROR FUNCTION* 

By Walter Gautschi^ 

Introduction . The error function of a complex variable, in more or 
less disguised forms, occurs in many branches of science and technology. 
Properties of this function, and computational methods, have been 
studied extensively. A useful survey, as of 1966, may be found in [1], 
and more recent work in [2], [11] , [14] . In many applications the func- 
tion must be evaluated a large number of times. It is therefore important 
to search for methods which are as efficient as possible. Current practice 
attempts to achieve the desired economy by adopting different methods in 
different regions of the complex plane. In this paper, instead, we pro- 
pose a single algorithm which is uniformly effective for all complex argu- 
ments. A corresponding ALGOL procedure is to appear in [8], 

In section 1 we review some relevant mathematical properties of the 
complex error function. Although much of this material is known, a few 
remarks are made which do not seem to be common knowledge. Among these is 
the observation that a certain continued fraction, known as the Laplace con- 
tinued fraction, while divergent on the real line, approximates the error 
function asymptotically in the sense of Poincare. The computational algo- 
rithm is developed in section 2. Basically, it consists of evaluating a 
truncated Taylor expansion. The increment, h, as well as the number of 
terms, N, are made to depend on the argument z at which the function is 
evaluated. As |z| increases, h and N decrease, until eventually both 

__________ 

Work supported, in part, by the National Aeronautics and Space 
Administration (NASA) under grant NGR 15-005-039 and, in part, 
by Argonne National Laboratory. 

^Department of Computer Sciences, Purdue University, Lafayette, 
Indiana 47907. 
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become zero, at which time the algorithm reduces to that of evaluating the 
Laplace continued fraction. Some performance characteristics, and data on 
testing the algorithm, are included in section 3. 


1. Mathematical preliminaries. The function 


_ 2 

( 1 . 1 ) w(z) = e Z erfc(-iz). 


where erfc ? 


(2/Vtt) 



dt denotes the complementary error function, 


was first introduced (and tabulated) by Faddeeva and Terent’ev [4], As a 
function of the complex variable z, w(z) represents an entire function, and 
has the property that both its real and imaginary parts are between zero and 
one for z in the first quadrant of the complex plane. This property may 
well have been one of the motivations for considering ( 1 . 1 ) as the basic 
form of the error function for complex argument. 

Closely related to (1.1) is the integral 


( 1 . 2 ) 


f(z) 



* —00 


Z-t 


dt. 


We have in fact 
(1.3) w (z) = 


i f (z) 

TT 

2 

7 f(z) + 2 e~ Z 

IT 


(Im z > 0) 
(Im z < 0) . 


While w(z) is an entire function, f(z) is analytic for all z not on the 
real line, and represents two analytic functions , one in the upper, another 
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in the lower half-plane, neither of which is the analytic continuation of 
the other. For real z, the integral in (1.2) is meaningful only in the 
sense of a Cauchy principal value integral. 

We note that (1.2) is a special case of a Stieltjes transform 

*00 

da(t)/(z-t). Most of the properties to be described below follow 

' —00 

from Stieltjes' classical theory [12], [10] concerning integrals of this 
type, and are therefore applicable in other situations as well (e.g», in 
the computation of the complex exponential integral) . 

Expanding the integrand in (1.2) in descending powers of z, and in- 
tegrating term by term, one obtains the asymptotic expansion 


(1.4) f(z) 'v l 


s=o z 


s+1 


(z » in |lm z| >_ a, a > 0 arbitrary), 


where 


(1.5> p s = 


• 2 
s -t 

t e dt = 


0 (s odd) , 

F((s+l)/2) (s even) 

2 2 
“”t —2 

are the moments of e . Since e has the zero asymptotic expansion in 

and SJl/4.4 < ffSr /4 

-tt/4 < arg z < fa/kj it is not surprising, in view of (1.3), that 


(1.6) w(z) ^ — V 

IT *•* 


s=o z 


s+1 


(z “ in -ir/4 < arg z < 5 tt/4) 


With the (formal) expansion (1,4) is associated the continued fraction 



z- 


1/2 1 _ 
z- z- 


3/2 

z- 


(1.7) 


a q o 


9 
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known as the Laplace continued fraction. More precisely, (1.7) is associated 
with (1.4) in the following sense. Let 


( 1 . 8 ) 


% Ill i. 1/2 (n-l)/2 q n (z) 

z- z- z- z- *** z p (z) 

n 


denote the n-th convergent of the continued fraction (1.7). It is easily 
verified that q n (z) is a polynomial of degree n-1, while p fl (z) is a monic 
polynomial of degree n. Then the quotient in (1.8) , if expanded in descend- 
ing powers of z , 


(1.9) 


vw 

Pn (z> 


I 


s=o 



s 


s+1 

z 


» 


yields a power series which agrees with that in (1*4) up to and including 
_2 n 

the term with z , i.e. 


(1.10) 


.(n) 


for s = 0,1,2, ... ,2n-l. 


It is also known [13] that the continued fraction (1.7) in fact converges to 
f(z) for every nonreal z. 

Another remarkable connection of the continued fraction (1.7) with the 
integral in (1.2) is obtained if the rational function (1.8) is decomposed 
into partial fractions. 


( 1 . 11 ) 


V z> 



n 


I - 

k=l z 



Expanding both sides of this equation in descending powers of z, and compar- 
ing coefficients of like powers, one finds that 
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v (n) - I x' n) [t' n> ] S 

S k ‘i k k 

In view of (1.10) it follows that 


(s — 0,1,2,...) . 


V " I A< n) [ t < n) ] S for s = 0,1,2,... , 2n-l , 

k-1 k k 

which, on account of (1.5) , implies that and t^ n ^ are the weights 

and nodes, respectively, of n-point Gauss-Hermite quadrature. 
Consequently, 


( 1 . 12 ) 


n A, 

lim V — 

n-t~ k=l z-J n) 
k 



z-t 


dt 


(Irn z ^ 0) , 


i.e. , Gauss-Hermite quadrature, applied to the integral in (1.2), con- 
verges for every z not on the real line. 

Using the well-known remainder term of Gauss-Hermite quadrature it 
also follows that 


2n /iF n l 

•as - a 

2 n 

It is interesting to observe that, although (1.12) does not converge 

if z = x is real, the Gauss-Hermite quadrature sum (1.11) for fixed n and 

-2n-l 

Z m X °° nevertheless approximates -iir w(x) to within an error of 0(x” ) . 

In fact, this is true as z + » In the sector -tt/ 4 < arg z < 5ir/4. In other 
words, the quadrature sums (1.11) , and thus t/ze convergents of the Laplace 
continued fraction (1.7) approximate -iiT w(z) asymptotically as z -> 00 in 
-ir/4 < arg z < 5 tt/ 4. This follows by combining (1.6) and (1.9)-(1.11) , 


(1.13) 


2n 


(n) 

~ v 2n " 




n 


dt - 1 

k=l 
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If this is compared with the asymptotic expansion (1.6), i.e. 


one notes that the leading term on the right in (1.14) is smaller than the 
corresponding term in (1.15) by a factor of 



. 2n-2 

(1.15) w(z) - ~ j 

*TT “ 


R n! . /— %,,-n , . \ 

'b /it n 2 (n -*• «). 

2 n r (n+%) 


This is why Gauss-Hermit e quadrature is so much more effective for computa- 
tion than straightforward asymptotic expansion. 
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There is yet another approach to the continued fraction in (1.7), which 
involves the repeated integrals of the complementary error function. Con- 
sider (see [6] for notations) 


_ 2 

(1.16) w r (z) - e Z i n erfc(-iz) 


(n — 0,1,2,...), 


W_ 1 (z) 



» 


so that in particular 


(1.17) w q (z) “ w(z) . 

OO 

If In z > 0, the sequence {w n (z)} n __^is known to be a "minimal" solution of 
the linear second-order difference equation 


(1.18) 


y n+l 


iz 

n+1 


'n 


1 

2 (n+1) y n-l 


0 


(n *» 0,1,2, ...) . 


(For terminology, and subsequent development, see [7].) For any integer 
N >_ 0, and v > N, define 


(1.19) 


t 








m 1/2 

r n-l -iz + (n+l)r 

n 

v = r ..v i 
n n-1 n-1 


(n — v , v— 1 , ... ,0) , 


(n = 0,1,2, . . . ,N) . 


We shall write r^ V ](z), v^(z) for r , ,v , if we wish to indicate their 
n-1 9 n n-1 9 n 9 

dependency on v and z. It can then be shown that 


( 1 . 20 ) 


lim v^ (z) 
n 

V->°° 


w (z) 
n 


(Im z > 0, n » -1,0, 1,2,...), 
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and consequently, 

lim r^|(z) = w^(z) /w n _^(z) (Im z > 0, n = 0,1,2,...). 

In particular, by (1.17), 

(1.21) w(z) = lim v^(z) = — - lim r^(z) (Im z > 0). 

v O /— v —X 

\)-X» y*if V“*°° 

To see the connection with the Laplace continued fraction, let n = 0 
in the second line of (1.19), and then in turn n = 0,l,2,...,v in the first 
line of (1.19). One obtains 


„[v],. _ 1 1 1/2 1 3/2 v/2 

V o W (-iz)+ T=iz)+ (~iz)+ (-iz)+ (-iz) 

3 i 1/2 1_ 3/2 v/2 

~ ir z- z- z- z- * * * z * 


where the second expression follows from the first by an obvious equivalence 
transformation. Comparison with (1.8) shows that 


( 1 . 22 ) 



q v+l 

P v+1 


(z) 

Tz) • 


Curiously, the function w n (z) defined in (1.16) is related to the n-th 
derivative of w(z) by 


(1.23) w (n) (z) = (2i) n n* w n (z) 


(n 0,1,2,...), 


a result apparently first observed in (5, p, 223] 
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2. Computational procedure . Our objective is to devise an efficient 
procedure for computing the function w(z) defined in (1,1) to a given number 
d of correct decimal digits after the decimal point, i.e., to within an 
(absolute) error of 10 We shall assume z to lie in the first quadrant 
of the complex plane. This is no restriction of generality, since 

_ 2 _ 

(2.1) w(-z) = 2 e -Z - w(z) , w(z) « w(-z) 


can be used to continue w into the remaining quadrants. 

As shown in (1.12), (1.14), Gauss-Hermite quadrature, or equivalently, 
the Laplace continued fraction (1.7), provides an effective means of com- 
puting w(z) for z e and |z| large. To obtain a more concrete idea as to 
the errors involved, we construct the altitude map of the meromorphic 
function 


e n (z) = w(z) 


^ n 

-± j 

TT i-* 


v (n) 

V k 


U k=l z-t< n) 
k 


i.e., the curves of constant modulus |e (z)'| « r. These may be obtained by 
numerical integration of the differential equations 


dx 

d<j> 


-Im 


e n <2) 


& - Re 

e (z) 
n 

L e » <2) J 

9 

d<j> R 

e’(z) 

n 


(z = x+iy) , 


subject to the initial conditions 


x(0) = 0, 


y(0) = n. 
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where n is the root of | (iy) | = r. Selected results are shown in Fig. 1, 
where n = 9, and r = 10 -< *, d «* 2(2)10. 

Given any d, it is obviously possible to construct a rectangular region 

(2.2) Rj 0 <_ x <_ x q , 0 < y < y 

outside of which 9-point Gauss-Hermite quadrature yields an accuracy of d 
decimal places. For d = 10, e.g., Fig. 1 suggests the choice x q = 5.33, 
y Q = 4.29. Larger values of n would result in a smaller rectangle R, whereas 
smaller values of n would require a larger rectangle R for the same accuracy. 
It is not possible, in general, to arrive at an optimal choice of n, as such 
a choice would depend on the relative frequency with which the procedure is 
used in various regions of the complex plane. The choice n = 9 appears to 
be a reasonable compromise, and we shall fix this value for what follows. 

To compute w(z) outside of R, we thus apply (1.19) with v = 8 and N = 0. 
In view of (1.22), (1.11), this is indeed equivalent to evaluating w(z) from 
the integral representation (1.2) , (1.3) by a 9-point Gauss-Hermite quadrature 
rule . 

It remains to consider the case where z e R. If y = Im z is relatively 
small, a common procedure consists of computing w(z) from a Taylor expansion 
about z q = Rez, or alternatively, to write 

-z 2 2i 

w(z) = e + — F(z) , F(z) .« e 

/ iT 

and to expand F(z) about z q = Rez. There are two disadvantages to this 
approach; 
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(1) it requires the computation of Dawson's integral, F(x), for 
x ■= z q . Although good rational approximations are available for F(x) 

(see, e.g., [3]), the necessity of computing F(x) is apt to increase 
both the length of the program, and the total machine time, for com- 
puting w(z) ; 

(2) the recursive computation of the expansion coefficients is 
subject to considerable loss of accuracy, particularly for large x > 0. 

Interestingly enough, both these defects are removed if one expands 
"downward" rather than "upward", i.e. , if one computes w(z) from the 
Taylor expansion 

(2.3) „(z) f y . ( . n) ( - lh )», 

n • 

n=o 

where h > 0 is suitably chosen* This approach has the further advantage 
of being related to the Laplace continued fraction approach, and in fact 
gives rise to an algorithm which generalizes algorithm (1.19) (used outside 
of R). 

We observe from (1.23) that (2.3) can be written in the form 

00 

(2.4) w(z) = l (2h) n w n (z+ih). 

n=o 

Approximating w fl by v^ V ^ [cf. (1.20)], and truncating the infinite series, 
we are led to define 

a^ V ^(z,h) = l (2h) n vt^ (z+ih) . 
n=o 


(2.5) 
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Letting t = (2h) n v^ (z+ih) , one obtains from (1.19) the following 
n U 

algorithm to compute (2.5), 


(2.6) < t = o = — 

N ° ° fit 


1/2 

°> r n-l ~ h-i. + (,+l)r„ <■> * '>•'>-1 0), 

2 


n 


r -v 


t = 2h r .. t . , a = a , + t (n = 1,2,...,N). 

n n-1 n-l n n-X n 


We note that for h = 0, N = 0, algorithm (2*6) reduces to algorithm (1»19)* 
Given e = *— 10""^, it is clearly possible to determine N and v (both 
depending, on z , h, and e) such that 


(2.7) 1 ( z >h) - w(z) 1 £ e. 


Since the series in (2.4) converges for every z and h, we can indeed find 
N such that | (z ,h) - w(z) | £ e/2, and with N so determined, find v > N 
such that |a^ V ^(z,h) - (z,h) | £ e/2. The triangular inequality then 

yields the desired result (2.7). 

Efficiency being one of our major concerns, we propose to 
(i) let h, N, and v depend on z in such a way that h = N = 0, v = 8 
for z outside of R, the rectangle introduced in (2.2); 

(ii) empirically determine the smallest integers N and v, subject to 
(i) and compatible with (2.7), for each z e R. 

The motivation behind the first of these objectives is to arrive at a 
e-ingle algorithm for computing w(z) in all of Q^, viz. algorithm (2.6), 
which, as was already observed, automatically reduces to the Laplace 
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continued fraction algorithm (1.19) when h = N = 0. The objective can be 

attained in many different ways. Exploratory computations led us to set up 

( 2 ) 

h, N, v tentatively in the form 


h = h s(z), N = {N +N.s(z)}, v = {v -Hj-sCz)} if z e R, 

O 0.1 O 1 


( 2 . 8 ) 


h = N = 0, v «* 8 if z e Q-N^R, 


where 


s(z) = (1 - 2_) -\JT- 

•'o ’ 


(— ) 2 
X 

o 


(z » x+iy) 


and h o? N q , N^, v q , are parameters which remain to be determined. Our 
second objective (ii) will serve to determine the last four of these param- 
eters, while the first, h , will be chosen so as to minimize machine time* 

’ ’o’ 

A basic aid in this parameter study is a gauging routine r, which does 
the followings given z = x+iy, h, and e, it returns nearly optimal values 
N = Np, v = Vp compatible with (2.7). The detailed steps involved in T are 
as follows: 

A A 

(1) Select N = N sufficiently large (say, N = 40). 

v 

(2) Determine v, the smallest integer v > N such that 

v' 

max„ kI V+1 °^ (z,h) - cJj V ^(z,h)| <_ e/100. The quantities o^ V ^(z,h), 
0<N<N 

A 

N = 0,1, ... ,N, are considered sufficiently accurate to represent 
true partial sums of the Taylor series (2.4). 


m 


{u} denotes the integer closest to u, 
tained in u+1/2. 


i.e., the largest integer con- 
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(3) Find Nj, as the smallest integer N <_ If - 3 satisfying 

r*! r 

| CT jIg_j. 3 C z » h) - J (z,h)| _< e. If there is no such integer N, 

A 

increase N by 10, and repeat steps (l)-(3). 

(4) Determine v r as the smallest integer v satisfying 

l~ Iv+ l ] ( 2 ,h) - (z,h) | <_ e. 


N 


A first application of the gauging routine T is made in determining 
the parameter h Q . The choice of h Q affects both the convergence of 
a^ V ^(z,h), as v 00 , and the convergence of the Taylor expansion (2.4). 

In fact, large values of h Q give rise to fast convergence of a^ v ^(z,h), 
but slow convergence in (2.4), while small values of h Q yield slow con- 
vergence of 0 ^ v ^(z,h) (particularly if y *» Im z is small), but fast con- 
vergence in (2.4). A good choice of h Q is therefore one which strikes a 
balance between these two opposing effects. In order to search for such 
a value, we let y = 0 (where these effects are most pronounced) and apply 
T with input parameters z = x, h = h Q s(x), e, for selected values of x 
and h Q . After each application of T we measure the machine time required 
to compute a^ V ^(x,h) by algorithm (2.6), where v = v^,, N = N^, are the in- 
tegers returned by T. With x q = 5.33, e = JglO -10 , the results are shown 
in Fig. 2, where machine time (in milliseconds) is plotted versus h Q for 
x - 0(1)5. It is seen from these graphs that a good choice of h Q is 

h = 1.6. 
o 

We next apply T to determine the parameters N q , N^, v q , in (2.8). 

A first pair of constraints is obtained by letting z = 0 and requiring that 


N 

o 


+ N„ 



v + 
o 


o 


9 


(2.9) 
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o o 

where Np , Vp are the results returned by T with input parameters z = 0, 
h = h ', e. A narrow band near the separation line y = y Q is then examined 
more carefully, since our preliminary computations indicated that Np and 
Vp approach limits larger than 0 and 8, respectively, as y + y Q . With Np 

A 

the largest Np, and Vp the largest Vp returned by T for input parameters 
z (near the line y = y Q ) , h = h s(z)., e, we let 

(2.10) N q = Np , V q = Vp, 

which , together with (2.9), determines the desired parameters uniquely . 

In our case of interest (x q =5.33, y Q =4.29, h Q =1.6, e = %L0 ^) , the 

O O A ^ 

results, are Np = 29, Vp = 30, Np = 6, Vp = 9, giving 

(2.11) h = 1.6 s (z) , N = (6+23 s(z)}, v = {9+21 s(z)} if z e R. 

Note that v > N for all z, as required in algorithm (2.6). 

It remains to examine whether this choice of parameters is indeed com- 
patible with (2.7). This is done by applying F with input parameters z = x+iy, 
h ■ -h s(z), e over a grid of points z e R, and by checking the inequalities 

{N q + N 1 s (z) } > Np(z) , {v Q + v 1 s(z)} >_ Vp(z), 

where Np(z) and Vp(z) are the output values of T at the point z. Using the 
grid x = 0(.5)5(.05)5.4, y = 0(.2)4(. 05)4.4, and e = %L0 10 , it was found 
that both inequalities are consistently satisfied. 
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With the values of h, N, v, defined in (2.8) , (2.11), the desired func- 
tion w(z) is thus approximated by o^(z,h) in (2.5), which in turn can be 
calculated by algorithm (2.6). The result is essentially the Algorithm in 
[8], except that in this procedure the sum Oj| V ^(z,h) is evaluated somewhat 

differently. Letting s = [v^ (z+ih) ] _1 J (2h) k vj v ^(z+ih), s M = 0, one 

n n k-n+1 fc W 

can write indeed 


r =0, 
v * 


•k = °> 


(2 - 12 > < Vi - P 


1/2 


iz+(n+l)r 


n 


s •, = r . [ (2h) n + s ] 
n-1 n-l l ' n J 


and then has 


(2.13) oJj v] (* f h) - / 


St 


Sn 


-1 


-1 


(if n < N) 


n — v ,v— 1, ... ,0 j 


(h > 0) . 


(h = 0) . 


The advantage of this algorithm over algorithm (2.6) is its increased 
speed on a digital computer [cf. section 3] together with the fact that no 
array of storage must be provided to hold the quantities r n _^( n ** 1»2,...,N) 


3. Performance characteristics and tests . We begin by comparing the 
computing times of the algorithm in [8] (referred to below as Algorithm...) 
with those of a similar algorithm, in which (2.12), (2.13) is replaced by 
(2.6). Both algorithms are compiled and executed on the CDC 6500 computer. 
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We select five layers S n ( n = 0,1,..., 4) in R, defined by 

S n ! 0 1 x 1 x 0 » n y 0 /5 1 y 1 (n+l)y o /5 (x q = 5.33, y Q » 4.29), 

and time the algorithms on each S n> The average time on S n is obtained 

by measuring the computing time of evaluating w(z) for 1000 values of z, 

distributed uniformly in S^, and by dividing the measured time by 1000, 

Both algorithms are timed similarly in the region outside of R (where com- 

(3) 

puting time is independent of z) . The results are shown in Table 1. 

It is seen that the second algorithm is slower than the first by a factor 
of 1.6-2. 2. 


z in 

Computing time (in mill is ec.) 

Algorithm. . . 

(2.6) replacing (2.12) , (2.13) 

S 0 

6.7 

14.5 

S 1 

6.0 

12.6 

S 2 

5.2 

10.7 

S 3 

4.4 

00 

0 

■^J 

S 4 

3.6 

6.8 


2.2 

3.6 


Table 1. Timing of Algorithm. .. and a related algorithm. 

_ 

Such timings are subject to slight variations, even on the same 
computer, due to such incidental factors as compiler, executive 
system, clock reading routine, etc. 
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For comparison we also timed the library subroutine for the exponential 
function e for selected values of x in the interval 1 x <_ 20. The time 
observed was rather consistently .315 milliseconds. The computation of w(z) 

(both real and imaginary part) by Algorithm. . . thus takes about as long as 
7-21 exponentiations, depending on. the location of the argument z. 

In order to gain further confidence in the accuracy claimed. Algorithm... 
is run for x = 0(. 02)5. 32(. 005)5. 35, 5.4(,2)6, y ® 0(. 02)4. 28(. 005)4. 31, 

4.4(.2)5, the results being compared with those obtained by the same algorithm, 
where h is replaced by 1.6, N by 33, and v by 36. (This combination of param- 
eter values yields 14 correct decimal digits for z near the origin.) The 
largest absolute deviation is found to be 5.18 x 10 ^ in the . real part, and 

4.91 x 10 - ^ in the imaginary part. For the maximum relative deviation the 

—8 —9 

figures are 1.18 x 10 and 6,64 x 10 , respectively. 

The region very close along the real axis is known to be a difficult 

f 

region in which to compute w(z) , particularly its real part. To see how well 
our algorithm does in this region, it is used to compute — Re w(xH-iy) for 

/iT 

x ■ 0(1)10 and y = s*10 -r (s = 1,2, 3, 5, 7; r ** 4,3,2), The results are com- 
pared with those in Hummer's table [9] . Remarkably enough, although some of 
the answers (for large x and small y) have order of magnitude 10 ^ , there is 
agreement to 8 significant digits (the precision in [9]) , excepting occasional 
end figure errors of 1-7 units. 

Acknowledgments . The author is indebted to Professor Henry C. Thacher, Jr., 
for stimulating conversations, and to Mr. Thomas J. Aird for writing the pro- 
gram to produce the altitude map of Fig. 1. 
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